Analysis for Growth and calcification response of the common collector urchin, Tripneustes gratilla, to projected climate change: effects of warming and acidification.

rm(list=ls(all=TRUE)) 
results <- read.csv("Data/Urchin_ClimateChange_Data.csv", header=T, na="NA")

library("lsmeans") #post-hoc tests
library("effects") #plot effects of modeling
library("lmtest") #linear mixed modeling
library("lme4") #linear mixed modeling
library("lmerTest") #calculate p-values in models
library("emmeans") #post-hoc tests
library("cowplot") #arrange plots
library("MuMIn")
library("car") #levenes tests, Anova, qqPlots, 
library("tidyverse")
library("stats")
library("plyr")  #splitting, applying, and combining data
library("dplyr")
library("onewaytests") #allows for Welch test with unequal variances (spine lengths)
library("stats")
library("effsize")
library("FSA")
library("rcompanion")
library("broom") #turns model output into data frames
library("plotrix") #Calculates standard error of the mean

Analysis approach:

Build and run model
Check for normality of residuals
Check for homogeniety of variance of residuals
Look at model summary
Run anova to check for significance of main effects
Run post-hoc to test hypotheses

Data Assembly and Results:

Diameter Means:

Assemble Data for initial collection of urchins from hatchery

##Initial Collection = Day -24

InitialDiameter <-
  #seperate out the initial sizes of urchins on day -24, when desired conditions were met.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "-24") %>% 
      #sizes on day -24  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


InitialMean <- mean(InitialDiameter$Diameter)
InitialSE <- std.error(InitialDiameter$Diameter)

Results of linear mixed model approach

InitialMod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=InitialDiameter)
anova(InitialMod, type = 2)

## Type II Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq   Mean Sq NumDF DenDF F value Pr(>F)
## Temperature    0.0005035 0.0005035     1    19  0.1144 0.7389
## pH             0.0102269 0.0102269     1    19  2.3232 0.1439
## Temperature:pH 0.0003840 0.0003840     1    19  0.0872 0.7709

Check assumptions

#check assumptions
# 1. Normality of residuals
qqPlot(residuals(InitialMod))

## [1] 27  4

shapiro.test(residuals(InitialMod)) #no pass

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(InitialMod)
## W = 0.52108, p-value = 4.878e-11

hist(residuals(InitialMod)) #eh

# 2. Equal variances
leveneTest(residuals(InitialMod)~InitialDiameter$Treatment) #No pass

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  3  2.9592 0.04307 *
##       42                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(fitted(InitialMod),resid(InitialMod,type="pearson"),col="blue")

Assemble Data

##After acclimating before ramp up = Day-13

AcclimatedDiameter <-
  #seperate out the initial sizes of urchins on day -24, when desired conditions were met.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "-13") %>% 
      #sizes on day -24  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


AcclimatedMean <- mean(AcclimatedDiameter$Diameter)
AcclimatedSE <- std.error(AcclimatedDiameter$Diameter)

Results

AcclimatedMod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=AcclimatedDiameter)
anova(AcclimatedMod, type = 2)

## Type II Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    0.015013 0.015013     1    19  0.2154 0.64782  
## pH             0.254996 0.254996     1    19  3.6588 0.07097 .
## Temperature:pH 0.177404 0.177404     1    19  2.5455 0.12711  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#check assumptions
# 1. Normality of residuals
qqPlot(residuals(AcclimatedMod)) #not super normal but... see below

## [1] 42  2

shapiro.test(residuals(AcclimatedMod)) #passes

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(AcclimatedMod)
## W = 0.9577, p-value = 0.09326

hist(residuals(AcclimatedMod)) #looks normalish

# 2. Equal variances
leveneTest(residuals(AcclimatedMod)~InitialDiameter$Treatment) #Passes

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.3961 0.2573
##       42

plot(fitted(AcclimatedMod),resid(AcclimatedMod,type="pearson"),col="blue")

Assemble Data

##After ramp up and desired conditions are reached to begin experiment = Day 1

Day1Diameter <-
  #seperate out the initial sizes of urchins on day 1, when desired conditions were met.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "1") %>% 
      #sizes on day 1  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


Day1Mean <- mean(Day1Diameter$Diameter)
Day1SE <- std.error(Day1Diameter$Diameter)

Results

Day1Mod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=Day1Diameter)
anova(Day1Mod, type = 2)

## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    0.04177 0.04177     1    19  0.1922 0.66606  
## pH             0.38094 0.38094     1    19  1.7525 0.20128  
## Temperature:pH 0.89969 0.89969     1    19  4.1389 0.05611 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#check assumptions
# 1. Normality of residuals
qqPlot(residuals(AcclimatedMod)) #not super normal but... see below

## [1] 42  2

shapiro.test(residuals(AcclimatedMod)) #passes

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(AcclimatedMod)
## W = 0.9577, p-value = 0.09326

hist(residuals(AcclimatedMod)) #looks normalish

# 2. Equal variances
leveneTest(residuals(AcclimatedMod)~InitialDiameter$Treatment) #Passes

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.3961 0.2573
##       42

plot(fitted(AcclimatedMod),resid(AcclimatedMod,type="pearson"),col="blue")

Growth

Subset data frame to obtain growth measurements.

### GROWTH MODEL 1: using day 1 as the first official day of experiment when desired environmental conditions were met (24 days after getting urchins. In that 24 days, they were acclimated and then conditions ramped up)


Initial <-
  #seperate out the initial sizes of urchins on day 1, when desired conditions were met to standardize growth to this measurement.
  results %>% 
    select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "1") %>% 
      #sizes on day 1  
    select("Temperature", "pH", "Diameter", "TankID")

  mean(Initial$Diameter)

## [1] 16.03

  std.error(Initial$Diameter)

## [1] 0.2550688

Final <- 
  #seperate out the final sizes of urchins on day 126
   results %>% 
    select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>% 
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column
    filter(Day == "126") %>% 
      #filter for sizes on day 126
    select("Temperature", "pH", "Diameter","TankID")


Growth <- 
  #create column that is growth
  bind_cols(Initial, Final) %>% 
    #binds initial (Diameter) and final (Diameter1) sizes to calculate growth
  mutate(Growth = ((Diameter1 - Diameter)/Diameter*100)) %>% 
    #add column for growth calculation
  select("Growth")


resultsGrow <-
  #table of initial and final size data
  results %>% 
    select("Day", "Treatment", "Temperature", "pH", "TankID", "Diam1", "Diam2") %>% 
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
    filter(Day == "126")
resultsGrow <- bind_cols(resultsGrow, Growth)
#combines table of initial and final with growth information

#Figure out means of each treatment group
resultsGrow %>% 
  dplyr::group_by(Temperature, pH) %>% 
  dplyr::summarise(mean = mean(Growth), s.e. = se(Growth))

## # A tibble: 4 x 4
## # Groups:   Temperature [2]
##   Temperature pH         mean  s.e.
##   <fct>       <fct>     <dbl> <dbl>
## 1 ambient     acidified  327. 16.3 
## 2 ambient     ambient    323. 11.0 
## 3 heated      acidified  352. 15.0 
## 4 heated      ambient    376.  7.15

GrowthMean <- mean(resultsGrow$Diameter)
#total grpwth mean pooled for all treatments - not that necessary
Day1SE <- std.error(Day1Diameter$Diameter)
#same as above for SE

Results

#LMM for growth:
modGrow <-
  resultsGrow %>% 
  lmer(Growth~ Temperature * pH + (1|TankID), data=.)

#Need to decide on model type that we want to use. 

summary(modGrow)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 417.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.90063 -0.32406  0.08141  0.32903  1.94882 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 1929.2   43.92   
##  Residual              283.1   16.83   
## Number of obs: 46, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                  326.595     18.578  19.000  17.580 3.28e-13
## Temperatureheated             25.439     26.273  19.000   0.968    0.345
## pHambient                     -3.424     26.273  19.000  -0.130    0.898
## Temperatureheated:pHambient   27.731     38.073  19.000   0.728    0.475
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                      
## Temperatureheated:pHambient    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.707              
## pHambient   -0.707  0.500       
## Tmprtrhtd:H  0.488 -0.690 -0.690

anova(modGrow, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.

## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    1169.41 1169.41     1    19  4.1304 0.05634 .
## pH               74.92   74.92     1    19  0.2646 0.61290  
## Temperature:pH  150.20  150.20     1    19  0.5305 0.47527  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of normality of residuals and homogeneity of variance. Assumptions are met.

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modGrow))

## [1] 43 20

shapiro.test(residuals(modGrow)) #passes

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modGrow)
## W = 0.95491, p-value = 0.07263

hist(residuals(modGrow)) #looks normal

# 2. Equal variances
leveneTest(residuals(modGrow)~resultsGrow$Treatment) #Passes

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  2.1906 0.1033
##       42

plot(fitted(modGrow),resid(modGrow,type="pearson"),col="blue")

Display exploratory interaction plot.

#Interaction plot shows slight interaction
interaction.plot(resultsGrow$Temperature,resultsGrow$pH,resultsGrow$Growth,
                 type = c("l","p", "b", "o","c"), legend = TRUE,
                 trace.label = deparse(substitute(pH)),
                 fixed = FALSE,
                 xlab = deparse(substitute(Temperature)),
                 ylab = deparse(substitute(Growth)),)

Initial as day -13

Model growth using the start as the day conditions began to ramp up.

### GROWTH MODEL 2: If run the same model using day -13 as the first technical day of the experiment. This is the day conditions began to ramp up, so urchins were experiencing gradual increase of stressors.


Initial1 <-
    #seperate out the initial sizes of urchins on day -14, when conditions began to ramp up.
  results %>% 
    select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>%
    filter(TankID != "8t" ) %>% 
        #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
      #gather Diam1 and Diam2 into single column
    filter(Day == "-13") %>%  
      #filter for day at -13 to use this as the intial size.
    select("Temperature", "pH", "Diameter", "TankID")

  mean(Initial1$Diameter)
  std.error(Initial1$Diameter)


Final1 <- 
      #seperate out the final sizes of urchins on day 126
   results %>% 
    select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>% 
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
      #gather Diam1 and Diam2 into single column
    filter(Day == "126") %>% 
      #filter for size on last day of experiment
    select("Temperature", "pH", "Diameter","TankID")


Growth1 <- 
  #create column that is growth
  bind_cols(Initial1, Final1) %>% 
    #binds initial (Diameter) and final (Diameter1) sizes to calculate growth
  mutate(Growth = ((Diameter1 - Diameter)/Diameter*100)) %>%
    #add column for growth calculation
  select("Growth")


resultsGrow1 <-
  #table of initial and final size data
  results %>% 
    select("Day", "Treatment", "Temperature", "pH", "TankID", "Diam1", "Diam2") %>% 
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
    filter(Day == "126")
resultsGrow1 <- bind_cols(resultsGrow1, Growth1)
  #combines table of initial and final with growth information


#Figure out means of each treatment group
resultsGrow1 %>% 
  dplyr::group_by(Temperature, pH) %>% 
  dplyr::summarise(mean = mean(Growth), s.e. = se(Growth))

Analyze with linear mixed model.

#LMM for growth:
modGrow1 <- 
  resultsGrow1 %>% 
  lmer(Growth~ Temperature * pH + (1|TankID), data=.)


summary(modGrow1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 451.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.97009 -0.36352 -0.06635  0.38074  2.05538 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 4120     64.19   
##  Residual              667     25.83   
## Number of obs: 46, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                  556.515     27.244  19.000  20.427 2.17e-14
## Temperatureheated             56.177     38.529  19.000   1.458    0.161
## pHambient                     -3.451     38.529  19.000  -0.090    0.930
## Temperatureheated:pHambient    9.115     55.833  19.000   0.163    0.872
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                      
## Temperatureheated:pHambient    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.707              
## pHambient   -0.707  0.500       
## Tmprtrhtd:H  0.488 -0.690 -0.690

anova(modGrow1, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.

## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    3141.63 3141.63     1    19  4.7100 0.04287 *
## pH                0.68    0.68     1    19  0.0010 0.97489  
## Temperature:pH   17.78   17.78     1    19  0.0267 0.87205  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions.

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modGrow1))

## [1]  2 27

shapiro.test(residuals(modGrow1)) #passes

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modGrow1)
## W = 0.95274, p-value = 0.05986

hist(residuals(modGrow1)) #looks normal

# 2. Equal variances
leveneTest(residuals(modGrow1)~resultsGrow1$Treatment) #doesn't pass...

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)  
## group  3  2.9964 0.0413 *
##       42                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(fitted(modGrow1),resid(modGrow1,type="pearson"),col="blue")

Display exploratory interaction plot.

#Interaction plot shows interaction
interaction.plot(resultsGrow1$Temperature,resultsGrow1$pH,resultsGrow1$Growth,
                 type = c("l","p", "b", "o","c"), legend = TRUE,
                 trace.label = deparse(substitute(pH)),
                 fixed = FALSE,
                 xlab = deparse(substitute(Temperature)),
                 ylab = deparse(substitute(Growth)),)

Calcification

Tips of spines

Analysis of calcification at the tip of the spine:

Assemble data.

### Model 1: calcification at the tips of spines

resultsTip <- 
  #create data using only the SEM images of the tips of spines.
  results %>% 
  select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>% 
    #selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
  drop_na(RatioSEM) %>% 
  filter(Chosen == "yes", PartOfSpine == "Tip")
    #filters only for tips and chosen side. No dusty images to be analyzed.

Analyse with linear mixed model:

#LMM for calcification at the tips of spines 
modTip <-
  resultsTip %>% 
  lmer(RatioSEM ~ Temperature * pH + (1|TankID), data = .)

## boundary (singular) fit: see ?isSingular

summary(modTip)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 64.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.76413 -0.74507 -0.02052  0.63045  2.07630 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 0.0000   0.0000  
##  Residual             0.1356   0.3682  
## Number of obs: 67, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                  1.69118    0.08930 63.00000  18.938   <2e-16
## Temperatureheated           -0.06662    0.12453 63.00000  -0.535   0.5945
## pHambient                   -0.07462    0.12453 63.00000  -0.599   0.5512
## Temperatureheated:pHambient  0.30557    0.18089 63.00000   1.689   0.0961
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                      
## Temperatureheated:pHambient .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.717              
## pHambient   -0.717  0.514       
## Tmprtrhtd:H  0.494 -0.688 -0.688
## convergence code: 0
## boundary (singular) fit: see ?isSingular

anova(modTip, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.

## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    0.10158 0.10158     1    63  0.7492 0.39000  
## pH             0.08185 0.08185     1    63  0.6038 0.44006  
## Temperature:pH 0.38685 0.38685     1    63  2.8534 0.09612 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions.

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modTip)) #Normal

## [1] 38  8

shapiro.test(residuals(modTip)) #Pass

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modTip)
## W = 0.97737, p-value = 0.2605

hist(residuals(modTip))

# 2. Equal variances
leveneTest(residuals(modTip)~resultsTip$Treatment) #Pass

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.4545  0.715
##       63

plot(fitted(modTip),resid(modTip,type="pearson"),col="blue")

Conduct post hoc analysis.

## POST HOC ANALYSIS 
emmTip <- emmeans(modTip, ~Temperature*pH, adjust = "tukey")  #Estimated marginal means (Least-squares means)
multcomp::cld(emmTip) #create a compact letter display of all pair-wise comparisons

##  Temperature pH        emmean     SE   df lower.CL upper.CL .group
##  ambient     ambient     1.62 0.0868 18.2     1.43     1.80  1    
##  heated      acidified   1.62 0.0868 18.2     1.44     1.81  1    
##  ambient     acidified   1.69 0.0900 18.2     1.50     1.88  1    
##  heated      ambient     1.86 0.0993 16.6     1.65     2.07  1    
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05

pwpp(emmTip) #Constructs a plot of P values associated with pairwise comparisons of estimated marginal means.

Display interaction plot for exploration.

#Interaction plot for cacification at the tip of the spines - interaction
interaction.plot(resultsTip$Temperature,resultsTip$pH,resultsTip$RatioSEM,
                 type = c("l","p", "b", "o","c"), legend = TRUE,
                 trace.label = deparse(substitute(pH)),
                 fixed = FALSE,
                 xlab = deparse(substitute(Temperature)),
                 ylab = deparse(substitute(Growth)),)

Base of spines

Assemble data.

### Model 2: calcification in the base of the spines

resultsBase <-
   #create data using only the SEM images of the base of spines.
  results %>% 
  select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>% 
    #selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
  drop_na(RatioSEM) %>% 
  filter(Chosen == "yes", PartOfSpine == "Base")
    #filters only for tips and chosen side. No dusty images to be analyzed.

Results

#LMM for calcification at the tips of spines 
modBase <- 
  resultsBase %>% 
  lmer(RatioSEM~ Temperature * pH + (1|TankID), data=.)

summary(modBase)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 98.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5535 -0.6315 -0.1529  0.6918  1.9937 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 0.1303   0.3610  
##  Residual             0.1597   0.3996  
## Number of obs: 68, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                   2.1440     0.1749 18.8077  12.257 2.06e-10
## Temperatureheated             0.3566     0.2487 19.1707   1.434   0.1677
## pHambient                     0.6524     0.2474 18.8077   2.637   0.0163
## Temperatureheated:pHambient  -0.2224     0.3594 18.9804  -0.619   0.5434
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                   *  
## Temperatureheated:pHambient    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.703              
## pHambient   -0.707  0.497       
## Tmprtrhtd:H  0.487 -0.692 -0.688

anova(modBase, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.

## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Temperature    0.30743 0.30743     1 18.992  1.9251 0.181362   
## pH             1.48873 1.48873     1 18.960  9.3223 0.006553 **
## Temperature:pH 0.06114 0.06114     1 18.980  0.3829 0.543424   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions.

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modBase)) #Normal

## [1] 33 37

shapiro.test(residuals(modBase)) # Pass

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modBase)
## W = 0.97304, p-value = 0.1472

hist(residuals(modBase))

# 2. Equal variances
leveneTest(residuals(modBase)~resultsBase$Treatment) #Pass

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   1.068  0.369
##       64

plot(fitted(modBase),resid(modBase,type="pearson"),col="blue")

Post hoc analysis

## POST HOC ANALYSIS 
emmBase <- emmeans(modBase, ~Temperature*pH, adjust = "tukey")
  #Estimated marginal means (Least-squares means)
multcomp::cld(emmBase)

##  Temperature pH        emmean    SE   df lower.CL upper.CL .group
##  ambient     acidified   2.14 0.175 18.8     1.78     2.51  1    
##  heated      acidified   2.50 0.177 19.5     2.13     2.87  12   
##  ambient     ambient     2.80 0.175 18.8     2.43     3.16  12   
##  heated      ambient     2.93 0.192 18.8     2.53     3.33   2   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05

  #create a compact letter display of all pair-wise comparisons
pwpp(emmBase)

  #Constructs a plot of P values associated with pairwise comparisons of estimated marginal means.

Display interaction plot.

#Interaction plot for calcification at the base of the spines - no interaction
interaction.plot(resultsBase$Temperature,resultsBase$pH,resultsBase$RatioSEM,
                 type = c("l","p", "b", "o","c"), legend = TRUE,
                 trace.label = deparse(substitute(pH)),
                 fixed = FALSE,
                 xlab = deparse(substitute(Temperature)),
                 ylab = deparse(substitute(RatioSEM)),)

Dropped Spines

Assemble data.

## Number of loose spines counted at the bottom of tanks at the end of the experiment. These were either broken or shed. Note: these spines were not collected, just counted, so can not determine whether they were broken or shed.

resultsDropped <-
  #create data using only the needed variables.
  results %>% 
  select ("Treatment", "Temperature", "pH", "SpineCount", "TankID") %>% 
  drop_na(SpineCount)

Analyse with linear mixed effect model.

##LM of dropped spines - can't make it LMM because error: number of levels of each grouping factor must be < number of observations - only have one count for each TankID.
modDropped <- 
  resultsDropped %>% 
  lm(SpineCount~ Temperature * pH, data=.)

summary(modDropped)

## 
## Call:
## lm(formula = SpineCount ~ Temperature * pH, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0000  -6.0833  -0.1667   0.4000  20.0000 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   21.167      3.640   5.814 1.33e-05 ***
## Temperatureheated             -1.167      5.148  -0.227  0.82315    
## pHambient                    -21.000      5.148  -4.079  0.00064 ***
## Temperatureheated:pHambient    1.600      7.461   0.214  0.83248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.917 on 19 degrees of freedom
## Multiple R-squared:  0.6088, Adjusted R-squared:  0.547 
## F-statistic: 9.855 on 3 and 19 DF,  p-value: 0.0003899

anova(modDropped) #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.

## Analysis of Variance Table
## 
## Response: SpineCount
##                Df  Sum Sq Mean Sq F value    Pr(>F)    
## Temperature     1    1.52    1.52  0.0192    0.8914    
## pH              1 2345.78 2345.78 29.4995 3.062e-05 ***
## Temperature:pH  1    3.66    3.66  0.0460    0.8325    
## Residuals      19 1510.87   79.52                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##ASSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modDropped)) #not normal

## 1458 1460 
##    2    4

shapiro.test(residuals(modDropped)) #faaaiiilll

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modDropped)
## W = 0.86706, p-value = 0.005636

hist(residuals(modDropped))

# 2. Equal variances
leveneTest(residuals(modDropped)~resultsDropped$Treatment) # barely pass

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  3  2.8487 0.06483 .
##       19                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(fitted(modDropped),resid(modDropped,type="pearson"),col="blue")

Analysis does not meet assumptions. Attempt a square root transformation.

##DOES NOT MEET ASSUMPTIONS so ATTEMPT AT Transformation:
resultsDropped$tdata<-(resultsDropped$SpineCount)^1/2
  # transformation
modDropped1 <- lm(tdata~ Temperature * pH, data=resultsDropped)
  # run lm model of transformed data

###ASSUMPTIONS for Transformation
# 1. Normality of residuals
qqPlot(residuals(modDropped1)) #not normal

## 1458 1460 
##    2    4

shapiro.test(residuals(modDropped1)) #fail

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modDropped1)
## W = 0.86706, p-value = 0.005636

hist(residuals(modDropped1))

# 2. equal variances
bartlett.test(residuals(modDropped1)~resultsDropped$Treatment) #fail

## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuals(modDropped1) by resultsDropped$Treatment
## Bartlett's K-squared = 43.282, df = 3, p-value = 2.144e-09

plot(fitted(modDropped1),resid(modDropped1,type="pearson"),col="blue")

summary(modDropped1)

## 
## Call:
## lm(formula = tdata ~ Temperature * pH, data = resultsDropped)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0000 -3.0417 -0.0833  0.2000 10.0000 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  10.5833     1.8202   5.814 1.33e-05 ***
## Temperatureheated            -0.5833     2.5742  -0.227  0.82315    
## pHambient                   -10.5000     2.5742  -4.079  0.00064 ***
## Temperatureheated:pHambient   0.8000     3.7304   0.214  0.83248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.459 on 19 degrees of freedom
## Multiple R-squared:  0.6088, Adjusted R-squared:  0.547 
## F-statistic: 9.855 on 3 and 19 DF,  p-value: 0.0003899

Anova(modDropped1, type=2) #test for significance of model

## Anova Table (Type II tests)
## 
## Response: tdata
##                Sum Sq Df F value    Pr(>F)    
## Temperature      0.23  1  0.0118    0.9146    
## pH             586.44  1 29.4995 3.062e-05 ***
## Temperature:pH   0.91  1  0.0460    0.8325    
## Residuals      377.72 19                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Transfomraiton was not successful, so we use a non parametric Kruskal Wallis Test.

### Transformations not work, so use non parametric:
kruskal.test(resultsDropped$SpineCount~resultsDropped$Treatment)# SIG*** (p=0.0005563)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  resultsDropped$SpineCount by resultsDropped$Treatment
## Kruskal-Wallis chi-squared = 17.505, df = 3, p-value = 0.0005563

Conduct Dunn Post Hoc Test and display interaction plot.

### POST HOC Pairwise ######
dunn <- 
  resultsDropped %>% 
  dunnTest(SpineCount ~ Treatment,
           method = "holm", kw = TRUE, data = .)
dunnph <- dunn$res
cldList(P.adj ~ Comparison, data = dunnph, threshold = 0.05)

##   Group Letter MonoLetter
## 1   AMB      a        a  
## 2    OA      b         b 
## 3     T     ac        a c
## 4  T/OA     bc         bc

#Interaction plot for dropped spines - no interaction
interaction.plot(resultsDropped$Temperature,resultsDropped$pH,resultsDropped$SpineCount,
                 type = c("l","p", "b", "o","c"), legend = TRUE,
                 trace.label = deparse(substitute(pH)),
                 fixed = FALSE,
                 xlab = deparse(substitute(Temperature)),
                 ylab = deparse(substitute(SpineCount)),)

Analysis for Growth and calcification response of the common collector urchin, Tripneustes gratilla, to projected climate change: effects of warming and acidification.

Emily E. Sesno, Keisha Bahr, Ariana S. Huffmyer, Judith D. Lemus

Analysis approach:

Data Assembly and Results:

Diameter Means:

Assemble Data for initial collection of urchins from hatchery

Results of linear mixed model approach

Check assumptions

Growth

Initial as day -13

Calcification

Tips of spines

Base of spines

Dropped Spines